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^ Abstract 

^-^ Mitosis in the early syncytial Drosophila embryo is highly correlated in space and time, as manifested 

" in mitotic wavefronts that propagate across the embryo. In this paper we investigate the idea that 

,—1 the embryo can be considered a mechanically-excitable medium, and that mitotic wavefronts can be 

1-^ understood as nonlinear wavefronts that propagate through this medium. We study the wavefronts via 

»p-( both image analysis of confocal microscopy videos and theoretical models. We find that the mitotic 

^ wavefront can be resolved into two distinct wavefronts in each cycle, corresponding to metaphase and 

* ^ anaphase, respectively. The two wavefronts have the same speed and are separated by a time interval that 

is independent of cycle, supporting the idea that they are two different markers for the same process. To 

understand the wavefronts theoretically we analyze wavefront propagation in excitable media. We study 

two classes of models, one with biochemical signaling and one with mechanical signaling. We find that the 

^Lj dependence of wavefront speed on cycle number is most naturally explained by mechanical signaling, and 

r^ that the entire process suggests a scenario in which biochemical and mechanical signaling are coupled. 
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^^ 1 Introduction 

ly^ The early embryos of many species, including Drosophila (l}[3], Xenopus (4||6^, Oryzias M, Fundulus [8], 

^SJ and zebrafish l9J[lO], exhibit metachronous mitosis, in which mitosis progresses as a wavefront through 

^D the embryo. Such wavefronts are reminiscent of biochemical wavefronts that are used to transmit signals 

'^T across many cells in other biological systems, such as wavefronts of the molecule cAMP that propagate 

^^ in a colony of Dictyostelium when it begins to aggregate to form a fruiting body 11 13 . Propagating 

^^ wavefronts, however, need not be purely biochemical in origin. The process of mitosis is a highly mechan- 

^ ■ ical one that involves significant changes in the volume occupied by chromatin 14 as well as separation 



. . of chromosomes 15 . This raises the question of whether mitotic wavefronts are purely biochemical 

^ phenomena or whether they might have a mechanical component as well. 

k> The nuclei of the Drosophila embryo are syncytial (i.e., they share the same cytoplasm and are not 

^ separated into individual cells by plasma membranes) during their first thirteen division cycles. The 

^ nuclei migrate to the egg's surface during the ninth cycle. There they divide five more times, until the 

fourteenth cycle, when cell membranes form and gastrulation begins [I]. Mitotic wavefronts are observed 
in cycles 9 through 13 |1 . In this period, chemical diffusion is unhindered by membrane barriers. It 
is known that calcium, a signal carrier that influences many local phenomena including mitosis [16f 18 



exhibits localized waves or spikes of concentration in the syncytial embryo [19y-|23] . The observation by 
Parry et al. |20 that there is a calcium wavefront in each cycle that travels across the embryo at the same 
speed as the mitotic wavefront, and occurs near the onset of anaphase, suggests that calcium signaling 
may be responsible for the mitotic wavefront. 

However, mitosis is also a mechanical phenomenon. In the syncytial embryo, nuclei are embedded in 



an elastic cytoskelcton, which contains both actin and microtubules 24-26 . Actin caps assemble around 
each of the nuclei at the end of interphase, and provide anchor points for the mitotic spindles that pull 
the two daughter nuclei apart 124-27 . Recent work shows that mechanical interactions are important 



for re-organization of the nuclei after mitosis 28 , and optical tweezer experiments show that nuclei are 
mechanically coupled f29l. However, little is known about how mechanical interactions affect collective 
phenomena such as mitotic wavefronts at the level of the entire embryo. 

In this paper we report the results of both our image analysis of wavefronts in early Drosophila em- 
bryos, and our theoretical studies of models of wavefront propagation. Using novel tracking techniques, we 
analyzed confocal microscopy videos taken of Drosophila embryos in which the nuclear DNA/chromosomes 
are visualized by labeling their histones with GFP. Our analysis yields the position, shape and dynamics 
of the DNA/chromosomes with high temporal and spatial resolution during cycles 9-14. We observe 
two distinct wavefronts in each cycle, one corresponding to the onset of metaphase (at which point the 
chromosomes condense in the nuclear midplane, known as the metaphasic plate - see figure |4] for an illus- 
tration of the different stages) and one corresponding to the onset of anaphase. To our knowledge, the 
metaphasic wavefront has not been observed before. We find that for a given cycle in a given embryo, the 
metaphasic and anaphasic wavefronts have the same speed, which remains constant as they propagate 
across the embryo. Moreover, for a given embryo, the time between the two wavefronts is always the 
same, consistent with the view that the onsets of metaphase and anaphase can be viewed as different 
markers of the overall process of mitosis. Both wavefronts slow down with cycle, and produce large-scale 
collective motion of the nuclei. 

We treat the embryo theoretically as an excitable medium, consisting of nuclei that can be triggered 
into initiating metaphase or anaphase, thereby locally exciting the medium and thus signaling their 
neighbors. We not only consider the well-known case of nonlinear wavefront propagation in a chemi- 



cally excitable medium 30 31 , but introduce a model for the early embryo as a mechanically excitable 
medium 32 , through which mitotic wavefronts can propagate via stress diffusion. Comparing the data 
with the results of these two models, we find that our observations are difficult to reconcile with a purely 
biochemical scenario. In such a scenario, the wavefront speed has a tendency to increase with nuclear 
density, and thus with cycle, contrary to our observations. The observations can, however, be explained 
quite naturally by a novel scenario in which nuclei not only respond to their mechanical environment, 
but also actively use it to signal each other. Our results suggest that mitotic wavefronts in syncytial 
Drosophila embryos may constitute one example of a previously unexplored form of mechanical signaling 



via nonlinear wavefronts that could also arise in very different biological contexts 32 33 



2 Results 

2.1 Image analysis results 

Nuclear cycle and shape 

An example image of detected nuclei in a Drosophila embryo is shown in figure [l^. In each cy- 
cle, as the nuclei progress from interphase through metaphase to anaphase, the detected shape of the 
DNA/chromosomes changes in a well-defined manner (figure fib) . Newly separated nuclei are small and 
spherical, and thus show up in our shape tracking as small circles. During interphase, the nuclear DNA 
grows in size over time as it is duplicated. At the onset of metaphase, the chromosomes condense in the 
midplane of the nucleus, and appear to elongate into an ellipse. The final step of mitosis, the onset of 
anaphase, corresponds to two detectable changes in the shape: a sudden shift of the orientation axis over 
a 7r/2 angle, and a change of aspect ratio. An example plot showing the ratio of the length of the two 
axes as a function of time during a cell cycle is given in figure IT];. 



Metaphasic and anaphasic wavefronts 

The onsets of metaphase and anaphase, as determined by the axes ratio (figurellli) are indicated by dotted 
blue hnes and dashed orange hnes, respectively. Evidently the onset of metaphase forms two wavefronts, 
one propagating from each pole (not necessarily starting at the same time); the displacements of the 
chromosomes follow these wavefronts. Similarly, the onset of anaphase forms two wavefronts, also causing 
chromosomal displacements. Mitotic waves were first observed by Foe and Alberts [Ij; with better time 
resolution, it is evident that instead of a single mitotic wavefront in each cycle, there are two distinct 
wavefronts corresponding to the onsets of metaphase and anaphase. There may be additional wavefronts 
that cannot be resolved from the axes ratio or chromosomal displacements. 

Effect of shape changes on nuclear positions 

The processes of metaphase and anaphase affect not only the shapes of the chromosomes, but also their 
positions. After each of the shape changes, the nuclei move collectively through the embryo, almost 
exclusively along the long axis (which we designate as the x-axis), resulting in a global 'breathing mode' 
of the entire embryo (see SI movie 1 |34 ). Remarkably, after an initial transition in which the nuclei 
re-organize after anaphase, the nuclei hardly move with respect to their nearest neighbors during this 
collective movement. Figure [TJs shows the average displacement Ax along the x-axis of a small set of 
nuclei. Figure [TJ^ shows the same motion for all nuclei, illustrating that the metaphasic and anaphasic 
wavefronts are reflected in corresponding wavefronts in the mechanical response. 

Wavefront speeds 

We quantify the wavefront speeds in figure [2] for two sets of movies, where the environmental conditions 
(in particular the temperature) are approximately the same for all movies in a given set, but different 
from one set to the other (the data of the two sets were taken several months apart). Figure l2k shows an 
example of a position vs. time plot of all condensation (blue diamonds) and anaphase (red pluses) events 
in a single cycle. The slope, corresponding to the wavefront speed, is clearly constant across the embryo. 
Figure [2J3 shows the ratio of the speeds of the anaphasic and metaphasic wavefronts, showing that for a 
given embryo and cycle, the two wavefronts move at the same speed. From embryo to embryo there are 
large variations in wavefront speed (figure Pp), but they all show the same exponentially decaying trend. 
This trend is illustrated in figure [2J3, where we plot the same data, normalized by the speed of the first 
wavefront, on a log-linear scale. Figure [2ji shows that the time interval that separates the metaphasic 
from the anaphasic wavefront is the same for all cycles for a given embryo, but is different for the two 
different sets of data. By looking at the point at which the nuclear envelope breaks down and reforms. 
Foe and Alberts fl] also found that the duration of the mitotic phase is constant through cycles 10, 11 and 
12 (3 minutes in their observations, comparable to our result), but was longer for cycle 13 (5 minutes). 
The re-formation of the nuclear envelope membrane may therefore take significantly longer in the last 
syncytial cycle, even though the actual wavefronts keep following the pattern of the earlier cycles. 

Cycle statistics 

The nuclei on the surface are separated by a well-defined distance a„ , which decreases with cycle number 
n. Because the number of nuclei doubles from one cycle to the next, it is not surprising that a„ decays 
exponentially, scaling like a„ ^ 2^^^'^~"°\ with n the cycle number and no the number of the first 
observed cycle. We consistently found a value oi ^ — 0.46 in our experiments (figure [2f and table fl]). 
The value of /? is slightly less than 1/2, presumably because the curved embryo is being projected onto 
a plane. We have also measured the duration of each cycle, i„, and found that it increases with cycle 
number n, with a weak exponential growth: i„ — ioe°'^^ "i where {to — 33s for set 1 and to — 25s for 
set 2, see table [T] and figure [7|. 
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Figure 1. Observation of wavefronts and mechanical response, a) Image of a Drosophila 
embryo during mitosis at the end of cycle 11, with the detected chromosomal contours overlaid. 
Anaphasic wavefronts (orange dashed curved lines), the long axis (green dashed straight line) and a 
typical slice perpendicular to the long axis (green parallel straight lines) are indicated, b) Sketch of the 
three main states in image analysis: interphase (circular contours), metaphase (compressed elliptical 
contours), and anaphase (highly extended elliptical contours, perpendicular to metaphase contour). See 
also figure l4] c) Ratio of the two elliptical axes of the detected shape of the nuclear DNA/chromosomes 
vs. time in cycle 11, averaged over an a;-slice (as shown in a); error bars indicate variation within the 
slice. The transitions between interphase and metaphase, as well as the onset of anaphase, are sharp and 
indicated respectively by dotted (blue) and dashed (orange) vertical lines. The slice shown was taken at 
X = 200/xm. d) Kymograph showing the elliptical axes ratio, a/b (where white indicates values larger 
than 1 and black indicates values smaller than 1), as a function of position x and time. The dotted and 
dashed lines indicate the two wavefronts. e) Average x-displacement Ax of the nuclei within one slice 
vs. time. After a nucleus has divided, we use the average position of its two daughters. The slice shown 
is identical to the one in figure c. f ) Kymograph showing the collective motion of nuclei in slices taken 
at different positions along the long axis of the embryo. White indicates motion in the positive x 
direction, black in the negative x direction. Dotted and dashed lines again indicate the two wavefronts. 



2.2 Theoretical Analysis 

Our observation that the metaphasic and anaphasic wavefronts propagate at constant speed across the 
embryo suggests that the embryo can be considered as an excitable medium that supports nonlinear front 
propagation. Alternatively, the nuclei could all have biological clocks that determine when mitosis starts, 
which operate independently; in that case the wavefront would be only a result of a lucky timing of those 
clocks. We briefly discuss various timing models below, and show that these are inconsistent with our 
observed data, so we conclude that the nuclei do communicate. We have considered two distinct classes of 
models for front propagation in excitable media. In the first model the nuclei communicate by releasing 
a small chemical, which then diffuses to neighboring nuclei, triggering them to initiate metaphase or 
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Figure 2. Wavefront propagation and speeds, a) x-coordinate of nuclei at the onset of metaphase 
(blue diamonds) and anaphase (red pluses) vs. time for the wavefront shown in figure [I] Both events 
show two clear wavefronts moving in from near the embryo poles (solid lines), b) Ratio of the speeds of 
the metaphasic and anaphasic wavefronts for different embryos and cycles. Each embryo is indicated by 
a different symbol and color, with the closed and open symbols representing two different measurement 
sets. Ratios for a given cycle and different embryos are slightly separated horizontally, c) Wavefront 
speed vs. cycle. Two of the embryos contribute two waves per cycle (coming in from opposite poles, as 
in figure ^; blue squares and green diamonds) . Although the actual propagation speeds vary 
significantly from one embryo to the next, they all follow the same trend, decreasing with successive 
cycles, d) Time interval between the metaphasic and anaphasic wavefronts vs. cycle, e) Log-linear plot 
of wavefront speeds vs. cycle, normalized by the speed of its first occurring wave front (if the first wave 
front is in cycle 10) or 0.71 times its first occurring wavefront (if the first wave front is in cycle 11). The 
black open circles connected by a dashed line corresponds to a scaling of 0.71 per cycle, showing that all 
embryos follow the same exponentially decaying trend, f) Average distance between nearest neighbors 
on a logarithmic plot. The dashed line has a slope of 0.46. In figures b-f, the same symbol/color 
corresponds to the same embryo. 

anaphase. In the second model we explore the novel idea that mitotic wavefronts in the early embryo can 
be described by wavefront propagation in a medium that is mechanically rather than chemically excitable. 
In this model, the deformation of nuclei during metaphase or anaphase gives rise to mechanical stresses 
that trigger other nuclei to proceed to that phase as well. 



Timing models 

Consider the following timing mechanism for generating a wavefront in a row of people. Assume each 
person has a (synchronized) watch, and each is told to raise his/her arms at a constant time interval apart. 
If the clock of each successive person in the row runs a little slower than that of the previous person, and 
they all raise their arms when their local clock reaches a certain agreed upon time, a wavefront will be 
generated. 



For the timing method to be the cause for the wavefronts observed in our system, each nucleus would 
require a clock. That clock could simply be the amount of time it takes to duplicate the DNA, i.e., the 
duration of the interphase, which changes from one cycle to the next. However, there is no correlation 
between a nucleus' position and the duration of interphase, which means that we would not expect a 
mitotic wavefront to emerge in this case. Alternatively, it is well-established that there are several proteins 
which exhibit patterning along the anterior-posterior or dorsal- ventral axes of the embryo. A much studied 
example is Bicoid, which exhibits an exponential profile along the anterior-posterior axis 35p6 , the same 



axis along which the wavefront travels. Now if the duration of interphase were affected by the local Bicoid 

concentration, that could provide a mechanism for the clocks of the nuclei to get out of sync, and produce 

a wavefront in the various markers for mitosis (such as our observed metaphase and anaphase wavefronts) . 

There are three reasons why the model outlined above cannot explain our data. The first is specific to 



Bicoid. As observed by Gregor et al. 36 , the total amount of Bicoid steadily increases over time, as more 
of the protein is translated in each cycle. In particular, the amount of protein keeps steady pace with 
the number of nuclei, such that at the start of each cycle, the actual amount of protein in each nucleus 
at a given position in the embryo is always the same. Therefore if Bicoid were responsible for causing 
the mitotic wavefronts, the wavefronts would have the same speed in each cycle, which they do not. The 
second reason is more general: as shown in appendix|D] in order to obtain a linear wavefront propagation 
from an exponential concentration profile, the actual absolute amount of material does not matter, only 
the decay length - which means that once again the predicted wavefront speed would be independent of 
cycle. Finally note that to get wavefronts traveling in both directions along the anterior-posterior axis of 
the embryo, we would need at least two concentration gradients of different proteins, for which it would 
be highly surprising if they produced wavefronts with the same speed. The timing method therefore 
cannot describe our data. 

Biochemical-signaling model 

At the end of a cycle, when all nuclei have completed the duplication of their DNA, we assume that 
they are in an excitable state, meaning that they can be triggered to initiate mitosis once they receive 
an appropriate signal. An obvious candidate for signaling between nuclei is the diffusion of a small ion, 
molecule or protein {e.g. calcium, cyclic AMP, Bicoid), which we will denote as A. By definition, nuclei 
can divide only once per cycle; therefore, in our model, we introduce a refractory period for each nucleus 
following anaphase, equal to the duration of the interphase. 

To introduce chemical excitability, we assume that if the local concentration of A exceeds a threshold 
a, the nucleus starts its program of mitosis, part of which involves releasing more A. A then diffuses 
away, raising the concentration of A at neighboring nuclei, and so on. In our model we also allow for 
a time delay idciay between trigger and release, meaning that a nucleus does not release more A until a 
time idciay after its local concentration exceeds a. We model releases of A by the nuclei (or sources) as 
localized pulses (Dirac delta functions), and the system is initiated with a single nucleus releasing A. The 
wavefront at any point in time corresponds to the position of all nuclei that release A at that moment. 
Details on how to solve the diffusion equation and carry out the other needed calculations are given in 
appendices [E| and [G] An example wavefront is shown in figure [3k. 

The speed ■;; of the resulting wavefront is determined by three parameters: the diffusion constant D, the 
nuclear spacing a and the concentration threshold a. Because our system is effectively two-dimensional. 



a has dimensions of 1/length . We obtained the value of a from direct measurements. Gregor et al. 37 
found from diffusion experiments in Drosophila that the diffusion constant of a molecule with hydro- 
dynamic radius R is well described by a modified Stokes-Einstein relation 38): D — k^T / {QiirjR) -\- b, 
where k^ is Boltzmann's constant, T the temperature, t] = 4.1 ± 0.4cP the effective viscosity of the 
syncytial Drosophila embryo, and b = 6.2 ± 1.0/im^/s is an experimentally determined constant. Using 
this expression, we estimate that a calcium ion with a radius of approximately 0.5nm has a diffusion 
constant of about 10^/im^/s. We have no direct way of estimating the threshold a, so we solved the 



model numerically for a range of values of a. 

Combining the parameters of our model, we define a nondimensional threshold, a = Aira'^a and speed 
V = -jy/-- As shown in appendix p^]| we then have v — l//(a), where /(a) increases monotonically 



with a (figure fl|. Consequently, TFboth D and a are fixed, the wavefront speed v increases as the 
nuclear spacing a decreases, and thus the speed increases with cell cycle, in direct contradiction to our 
experimental observations. Moreover, numerically solving for the wavefront in this model shows that 
the resulting wave speeds are much larger than the ones we observed in the experiments, by up to two 



orders of magnitude in the later cycles (see figure 10). For a different chemical with a smaller diffusion 



constant we find that the magnitude of the predicted speed is closer to the experimental value, but the 
trend (increasing speed with cycle) persists. Thus, the simplest form of the biochemical signaling model 
cannot describe the data of figure [2J:. 

We next consider the possibility of a delay idciay between the time when the local concentration of A 
reaches the threshold value a, and the instant when more A is released. For details on the implementation 



see appendix G.2 In the limit where c? jD ^ idciay, the wavefront speed is determined by diffusion as 
before, v = D/{af{a)). In the opposite limit, a^/D ^ idciay, we find v = a/tdciay, so v would decrease 
with cycle number for constant tdciay However, we find that to get the correct speeds in the earlier cycles 
requires a delay time of 5-lOs, which is comparable to a^ /D for calcium; for this range of delay times, the 
speed still increases with cycle number (blue diamonds in figure pp). Only for much larger delay times 
of 100s or more do we arrive at the regime where v — ajt^^x^y^ but then the speeds we obtain are much 
lower than those we measured experimentally. Again, for a chemical with a smaller diffusion constant, 
the trend remains the same; in that case o? / D is larger, and we move further into the diffusion-dominated 
regime. Thus, a biochemical-signaling model with a time delay that is independent of cell cycle cannot 
describe our observations either. 

We also investigated the wavefront speed in the case where the delay time is allowed to vary from 
one cycle to the next. Since the cycle duration shows a weak exponential increase with cycle number 
(figur e [7| , we might speculate that the delay time also increases exponentially. However, as shown in 
figure [To] (red points), this dependence cannot explain the experimental data either. As stated above, 
a diffusion constant D = 10^/im^/s requires a delay time of 5-lOs to fit the wavefront speed in the 
first observed cycle; however, the last cycle then requires a delay time of approximately lO'^s, which is 
comparable in length to the entire cycle and much longer than the duration of mitosis. We therefore 
consider smaller diffusion constants, although those would imply larger signaling molecules which have 
not been identified. In the case that D = lO/xm^/s (almost the lower limit for which we can still fit the 
speed of cycle 10), we can fit the data only if we allow the delay time to vary freely from one cycle to 
the next (see table l2|. However, we do not find a systematic dependence of the delay time on the cycle; 
moreover, even though the required delay time is now shorter than the duration of mitosis, in the latter 
cycles it is still significantly larger than the time it takes the wavefront to travel from one nucleus to the 
next. Consequently, whereas in the early cycles we would have simple neighbor-neighbor interactions, 
in the latter cycles there would be interactions between nuclei that are 4 or 5 rows apart. More details 
on this are given in appendix |E] There, we also analyze the possibility that the concentration threshold 
changes with cycle, and find that in order to fit the data, it has to change by four orders of magnitude. 
Even changing both idciay and a does not resolve these problems (tableful). 

On the basis of these results, we conclude that a wavefront that propagates via diffusion of calcium 
should not slow down with cycle number. It is possible that a wavefront of a larger biomolecule with 
a lower diffusion coefficient could do so, although this is also unlikely and no such species has been 
identified. We also note that any model in which the biochemical signal is mediated by a method that is 
faster than diffusion (such as active transport) suffers from the same problem: the predicted wavespeed 
would go up with increasing cycle, because the spacing between the nuclei goes down. 



Mechanical-signaling model 

The early embryo cannot support ordinary elastic waves because it is heavily damped by the viscosity 
of the cytosol. Consequently, displacements do not propagate ballistically as in a wave, but diffusively. 
However, just as diffusion of A can lead to wavefront propagation in the biochemical signaling model, 
diffusion of displacement could lead to wavefront propagation in a mechanical signaling model. We 
therefore introduce a model in which the nuclei communicate via stresses or strains that they exert on 
the cytoskeleton when they undergo condensation or anaphase. In this model, a nucleus starts its mitotic 
program when the largest eigenvalue of the local stress tensor exceeds a threshold value a. We describe 
the cytoskeleton as a homogeneous linear elastic medium, characterized by two elastic parameters, for 
example its bulk and shear moduli {K and /z, respectively) or the Young's modulus E and dimensionless 
Poisson ratio v. The viscous fluid in which the elastic cytoskeleton is immersed exerts a drag force on it, 
characterized by a damping constant T. Assuming that the nuclei exist in a thin layer near the surface of 
the embryo, we we denote the deformations in the plane of the layer by Ui — x'^^ — Xi (i = 1, 2), when the 
deformation maps point (xi,X2) onto point (x']^,a;2). The displacement u of a nucleus can be described 



by 39 



IT' /T" 

The term on the left represents the damping with damping factor F, and the two terms on the right are 
the elastic force. Equation (fTl) is reminiscent of the diffusion equation: a time derivative on the left equals 
second-order space derivatives plus a source term on the right. This model can therefore be thought of 
as describing the diffusion of the vector displacement field Ui . The right hand side of equation (IT]) gives 
rise to two quantities with the dimensions of diffusion constants f32' : 

„ E 1 — V li ,„ E LL 

Di = - ^— = and D, = — ^ = -• 2) 

In order to introduce mechanical excitability into the model, we assume that if the largest eigenvalue 
of the stress tensor at a nucleus at position xq exceeds a threshold value, a, at time to, it triggers a force 
dipole, adding another term djQijS{x — xo)Q{t — to) to the right hand side of equation (nj): 

E E 

TdtUi = ^7T~^^i^o'^' + , _ . didjUj + djQijS{x - xo)id{t - to), (3) 

where Qij is a symmetric tensor of rank 2. Here 6{x) is the two-dimensional Dirac delta function and 9(t) 
is the Heaviside step function. Equation ([s]) essentially describes the diffusion of the vector displacement 
field Ui due to a tensor source term. It is similar, but not identical, to a scalar reaction-diffusion equation, 
which describes the evolution of a scalar concentration field c due to a scalar source term. It is therefore 
not surprising that the model described by equation (fTl) also produces wavefronts, as can be seen in 
figure [3]3. 

In order to compare the model results with the data, we need to estimate the values of the elastic con- 
stants and the damping parameter. The speed v now depends on the quantity D = fj,/T that determines 
the dimensional part of both diffusion constants (equation ([2])), as well as the nuclear spacing a and the 
threshold value a. Following Palmer et al. pO], we estimate the shear modulus of our elastic network to 
be of the order of I — lOPa. The damping parameter F is the drag per unit volume, which for a mesh of 
semiflexible polymers should be roughly equal to the number of filaments times the drag on each of them: 
F ~ cryZmcsh '^ '7/^mosh' wherc c is the filament concentration, /mesh the mesh size and r] the viscosity of 
the surrounding fluid. For ry ~ 4 • 10~^Pa • s and ^mesh ^ lOOnm, we find that D '^ 5 — 20/xm^/s. 

The elastic moduli and the damping parameter should depend both on frequency (which is low in 



our case) and local filament concentration 40-44 , which can differ from one cycle to the next. Because 



the number of nuclei doubles in each cycle, the number of actin caps making up the network doubles as 



well (see figures [4] and Isl) . Thus, the local concentration of actin and of microtubules should effectively 
double with cycle number n. We therefore write c ~ 2'"~"°\ where as before rig is the number of the 
first observed cycle. Both the storage and loss moduli of polymer networks increase with concentration 



approximately as power laws, but the actual powers are debated [40[|42 44 . Moreover, in each successive 
cycle the nuclei get pushed further out into the plasma membrane encompassing the entire embryo \U, 
increasing the friction coefficient. Because the dynamics of our system depends only on the value of the 
two effective diffusion constants given in equation (pi) , we will not be able to distinguish the dependence of 
the storage and loss moduli independently. Instead we assume a dependence D = n/T ~ c~'^ ~ 2~t("~"'o). 

Because of the mathematical similarity between the mechanical-signaling model (equation (II])) and 
the diffusion model for concentration fields, we use the same type of dimensional analysis as for the 
diffusion-signaling model and write v — [D /a)g{a,v) — {^/aT)g{a,v), where now a = a?a/Q, with Q 
the magnitude of the force dipole and a the spacing between the nuclei. We determine g(a, v) numerically, 
and find that it can be well described by the functional form g{a, v) — —4(1 — a) log(a)/(l — v) [32]; for 
details, see the appendix |G.3| We now consider two possibilities. First, the value of the dimensionless 
threshold a can be constant throughout stages 9-14, which means that the individual nuclei are effectively 
the same in all cycles, adapting for the change in polymer density and available area. It is well known in 
the literature that such adaptive patterns are present in, for instance, expression patterns of morphogen 
proteins like Bicoid [36] , where the total amount of Bicoid increases in such a manner that the amount per 
nucleus remains the same. It is conceivable that something similar happens with the stress threshold a, 
causing the dimensionless version a to be constant. If this is the case, the wavefront speed only scales 
with the prefactor fi/oT, which itself depends on the cycle as 2^J^~'^'>^^~^°\ using the experimental result 
that a = Co • 2^'^'^"^"°). For this case we find from fitting the experimental data that 7 — 1.0 (green 
fits in figure prl). The second possibility is that the actual stress threshold value a is constant, and thus 
a scales like 2~2'^("~"<'). We then use our functional form for g{a,i') to fit the experimental data, and 
find that 7 = 1.4 and a — 0.1 • Q, where Q is the strength of the force dipoles (blue lines in figure l3]i). 
Note that in both cases we need to set the absolute scale, i.e., the wavefront speed in the first fitted 
cycle; this value vq absorbs the values of /^, F and a in the first fitted cycle as well as the factor 1 — ;/ in 
the denominator of our functional form. However, once vq is determined for each dataset, the same fit 
parameters (7 for fixed a and 7 and a for fixed a) describe both sets of data, indicating that the trend 
is the same. In both cases (a constant and a constant) we find a good fit; we can therefore at present 
not distinguish between the two possibilities. Further experiments, such as ones in which the embryo is 
poked mechanically with a well-controlled force (see Discussion), are needed to resolve this issue. 

The value for 7 obtained above is reasonable for crosslinked semiflexible networks. For a semidilute 



solution of rigid rods, the viscosity is expected to rise as C^, where c^ is the filament concentration 45 

j-2 j-3 "— 

'mosh'ent 



The modulus of a semiflexible network should scale as ^„„„i,^ent; where /mesh is the mesh size and Zent is the 



filament entanglement length 461. For a crosslink concentration c^^ , one would expect Imosh ~ ^cnt ~ Ca; , 
leading to /z ~ Cx ■ For the damping coefficient F we obtain F ~ CA^ll^lcsh ^ c\cx . In each cycle, 
we would expect Cx on ca, that is, the local concentrations of all biomolecules, actin and crosslinker 
proteins included, should increase proportionately. This would lead to -D = /Lt/F ~ c^ , so that 7 = 2, 
in reasonable agreement with 7 ~ 1.4. We note that it would be valuable to measure experimentally the 
actin concentration near the embryo surface as a function of cell cycle to pin down these estimates more 
precisely. 

3 Discussion 

During the early cycles of Drosophila development, the cycles of the nuclei are strongly coupled across 
the entire embryo by wavefronts of metaphase and anaphase onset that travel at constant speed across 
the embryo. We summarize our observations as follows: 
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Figure 3. Two models for the propagation of wavefronts by chemical and mechanical 
signaling, a) Color plot showing the chemical wavefront in two dimensions. The wave starts in the 
center (red dot) with a single Dirac delta peak release. The color coding indicates when a nucleus 
releases its chemical to the bulk, going from red through the different hues of the rainbow to violet, b) 
Color plot showing the mechanical wavefront in two dimensions, including the orientations of the 
dipoles, which are picked at random. The color coding is the same as in figure a. c) Plot comparing the 
different results of the two models to the two sets of experimental data (black and gray dots with error 
bars). Blue diamonds: biochemical model with fixed time delay of 6 seconds. Red diamonds: 
biochemical model with fixed time delay of 10 seconds. Green squares: mechanical model with 
D = /i/r ~ 2~'''("~"o) for 7 = 1.4 and a = 0.1 • Q. d) Fits of the scaling forms of the mechanical model 
(green and blue curves) to the experimental data (black and gray dots with error bars). Green lines: fit 
to w = Wo • 2'^''~'''^("~"''^ for the case where a is constant. Fit parameter 7 = 1.0. Blue lines: fit to 



v^vq- 2('3-'^)("-"«)g(a • 22/5("-»o)^ J,) where g{a, v) = -4(1 
threshold a is fixed. Fit parameters: a = O.IQ, 7 = 1.4. 



a) log(a)/(l — v), and the absolute 



1. There are at least two wavefronts in each cycle, corresponding to metaphase and anaphase, which 
travel at the same speed (figure l2b). 

2. The common speed of the metaphasic and anaphasic wavefronts slows down exponentially in each 
successive cycle (figures [2b andllK). 

3. The metaphasic and anaphasic wavefronts are separated by the same time interval in each cycle for 
a given embryo (figure pj: 



In addition to these observations, we add those of Parry et al. 20 



4. There is a calcium wavefront that coincides with the onset of anaphase. 

5. The speed of the calcium wavefront slows down in each successive cycle, matching the speed of the 
anaphasic wavefront. 

We now consider some scenarios to assess whether they are consistent with these observations. We 
note at the outset that the fact that the metaphasic and anaphasic wavefronts propagate with the same 
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speed, separated by a constant time interval, implies that there is a biochemical clock that sets the time 
interval. 

Scenario A 

A chemical wavefront triggers the metaphasic, calcium and anaphasic wavefronts. Here, we assume that 
there are clocks that cause both metaphase and anaphase to occur at specific time intervals after the 
biochemical signal triggers the clocks. In this case (4) implies that calcium cannot be the signal carrier, 
leaving the question what the signaling chemical is. Moreover, in this case, the results of the biochemical- 
signaling model suggest that (2) is very unlikely, since the model predicts that the speed of the wavefront 
should go up with cycle. Thus, Scenario A appears unlikely. 

Scenario B 

A mechanical wavefront triggers the metaphasic, calcium and anaphasic wavefronts. In this scenario, 
there is a mechanical wavefront that precedes the metaphasic wavefront. The stress released to maintain 
the mechanical wavefront starts a clock which triggers the onset of metaphase, calcium release and the 
onset of anaphase. Because the wavefront is mechanical, the wavefront speed slows down with cycle, 
according to the predictions of our mechanical signaling model. However, it is unclear what stress release 
triggers this wavefront, and how stress is released to amplify the signal to maintain the wavefront speed. 

Scenario C 

The metaphasic wavefront mechanically triggers the calcium and anaphasic wavefronts. In this scenario, 
the onset of metaphase itself at a nucleus leads to stress release. The stress starts a clock which triggers 
the release of calcium and also triggers the onset of anaphase. As a result, all three wavefronts travel 
at the same speed; because the metaphasic wavefront propagates mechanically, this speed slows down 
with successive cycles. In addition, the clock leads to fixed time intervals between the three wavefronts. 
This scenario is therefore consistent with all five observations. The scenario is also consistent with the 
expectation that the entire process of mitosis can be described by an extended wavefront, with metaphase, 
calcium release and anaphase as markers of this process that then also manifest as wavefronts. 

Scenario C is consistent as well with independent observations made in Xenopus embryos. These 
embryos are not syncytial; instead they are divided into cells from the first cycle. It is unlikely that a 
biochemical signal could cross cell membranes to propagate a wavefront. Nevertheless, these embryos do 
exhibit metachronous mitosis 4|. They also exhibit calcium oscillations inside each cell, which precede 



anaphase 47 . Their behavior is therefore most consistent with Scenario C: the metaphasic mechanical 
wavefront triggers a calcium signal inside each cell and an anaphasic wavefront. 

Scenario C implies several predictions. First, the metaphasic wavefront should precede the calcium 
wavefront, which precedes the anaphasic wavefront. Second, the speed of the calcium wavefront should 
be the same as the speed of the metaphasic and anaphasic wavefronts. Third, the time interval between 
the metaphasic and calcium wavefronts should be the same for all cycles for a given embryo. All of these 
predictions are consistent with the observations reported by Parry et al. [20), and can be further tested 
experimentally by combining labeled histones, as in our experiment, and calcium green dextran, as in the 
experiment of Parry et al. [20j . 

In principle, the estimated elastic constants and damping coefficients could be obtained directly from 
experiments by measuring the storage and loss moduli of the embryo surface in vivo using two-point 



microrheology. Optical tweezer experiments similar to the ones done by Schotz et al. 29 could also be 
used to extract the elastic moduli and the drag coefficient we used in our mechanical model. The actin 
concentration could be measured at the same time by staining the actin filaments with e.g. rhodamine. 



as done by Parry et al. 20 or GFP-moesin, as done by Cao et al. 48 
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The key assumption in Scenario C is that the initial metaphasic wavefront also propagates mechani- 
cally. This can be tested by mechanically poking the embryo at different times within the cell cycle. If 
the cell is poked just in advance of the metaphasic wavefront, the poking itself should generate a wave 
that propagates from the poking site with the same speed as the metaphasic wavefront. In addition, this 
poking wavefront should stimulate the release of calcium. One would therefore expect that metaphasic, 
calcium and anaphasic wavefronts could be generated by poking the embryo mechanically. If the embryo 
is poked too far in advance of the metaphasic wave, there should be no response because the DNA in the 
nuclei has not yet been fully duplicated, so mitosis is not yet possible. If the embryo is poked after the 
metaphasic wavefront begins, there may be no response because calcium would already be depleted by the 
metaphasic wavefront. Thus, we would expect that poking could generate an anaphasic wavefront only 
if it is applied in a certain time window of the cycle. Another method to test the mechanical-signaling 
model is to destroy the filaments that mechanically couple the nuclei, which should prevent the mechan- 
ical wavefronts from propagating and thus the nuclei from synchronizing their mitosis. This could be 
done by injecting colcemid to disrupt the microtubules or latrunculin which affects actin filaments [3]. 
Other means of disrupting cytoskeletal filaments, via mutation or laser ablation, should also affect the 
mechanical wave. 

4 Materials and Methods 

4.1 Confocal videos 

The imaged flies were from a His-GFP stock with a P[w+ ubi-II2A-GFP] insertion on the third chromo- 
some. All embryos were collected at 25° C and dechorionated in 100% bleach for 1 minute. They were 
picked using a 70fim nylon strainer (BD Falcon), rinsed in distilled water and laid down on a semiperme- 
able membrane (Biofolie) . The excess water was absorbed and the embryos were immersed in Halocarbon 
oil 27 (Sigma Aldrich) and covered with a 22 x 22^m coverslip (Corning). Embryos were imaged with 
a 20x oil immersion objective plan apochromat (Leica, NA=0.7) on a Leica SP5 laser scanning micro- 
scope with excitation wavelength of 488nm (argon laser 60mW) . 8 bit images were taken every second 
at 512 X 1024 0.45nm pixels and 1.4/is/pixel (734ms/image). An example video is shown in SI movie 1, 
available online [34) . 

4.2 Image analysis 

We visualized nuclear DNA/chromosomes by tagging their histones with GFP. To determine the positions, 
sizes, aspect ratios and orientations of the DNA/chromosomes from each video frame, we developed a 
new image analysis technique, explained in detail in [49, . In brief, we first applied a bandpass filter to 
eliminate high-frequency noise. We then made a contour plot of the resulting image, found the locally 
highest-level contour (i.e., the contours with no other contour inside them), and identified each of them 
as a single nucleus. For each nucleus, we fit the contour at half-height with an ellipse to get its position, 
shape and orientation. An example of an experimental image with the chromosomal tracking overlaid is 
given in figure [ik. 

Because the images were taken at high frequency (typically 1 Hz), the nuclei move less than their 
own radius from one frame to the next, simplifying tracking. The obvious exception is when nuclei divide 
during anaphase, and the observed shape splits in two. Because we detect shapes as well as positions of the 
chromosomes in each nucleus, tracking divisions is easy as well: when a nucleus divides, the chromosomes 
become highly elongated just before they split, and produce two almost circular daughters close to the 
endpoints of the long axis of the mother immediately after it splits, which are easily identified. 
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Appendices 

A Embryo layout and replication cycle 

The four main stages in the Drosophila embryo replication cycle, which we can detect from our movies, 
are illustrated in figure |4] A sketch of a cross-section of the embryo is shown in figure [5] illustrating how 
the nuclei are all located at the surface of the embryo for cycles 9-13 fTllSO]. 



Drosophila replication cycle 



interphase 



metaphase 



anaphase 



telophase 



interphase 




Figure 4. Illustration showing the four stages of the Drosophila embryo replication cycle that we can 
detect from our movies: interphase (DNA replication), metaphase (condensation of chromosomes in the 
nuclear midplane) , anaphase (division of the nucleus in two daughter nuclei) and telophase (separation 
of daughter nuclei) . The plasma membrane is shown in gray, the actin cap (made of actin filaments) in 
red, the microtubules in green, the centrosomes in yellow, and the DNA/chromosomes in blue. 
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Figure 5. Sketch of a cross-section through a Drosophila embryo valid for stages 9-13. Most nuclei are 
located at the surface of the embryo. The nuclei are pushed outwards into the plasma membrane (gray), 
resulting in the formation of somatic buds. Each nucleus is enclosed in a microtubule basket (green) 
and contained in an individual actin cap (red), which gets disassembled after mitosis and re-assembled 
during interphase. DNA/chromosomes are shown in blue and centrosomes in yellow. The yolk (light 
blue) is a viscoelastic fluid containing water, cytoskeletal elements and necessary building blocks for the 
nuclei. The yolk is bounded by an actin cortex over which the nuclei can move. Also shown in this 
sketch are the small number of nuclei that reside inside the yolk, and the also small number of somatic 
cells that already form in cycle 10 at the posterior end (the pole cells that divide out of sync with the 
rest of the embryo). See Foe and Alberts [l] for sketches for each of the first 14 cycles and Schejter and 
Wieschaus (50] for a review on the cytoskeletal elements in the early embryo. 
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B Experimental data sets 

Our image analysis results are for two different sets of experiments, which were carried out at ambient 
room temperature several months apart. The ambient temperature was higher for the second set, resulting 
in faster embryo development. We only used the data from those embryos which we could track from 
cycle 10-14 in the first set (Dataset 1, 3 embryos) and cycle 11-14 in the second set (Dataset 2, 4 
embryos). SI movie 1 (available online [34]) is the raw data of one of the embryos from set 1. This 
confocal microscopy imaging movie shows a developing Drosophila embryo. The chromosomal histones 
are visualized by labeling with GFP. The version of the movie shown here shows 1 image per 15s, displayed 
at 5fps, so sped up 75x. Movies for data analysis were recorded at Ifps. The dimensions of each frame 
are 346 x 440^m. 

C Additional image analysis results 

The average data from the two sets are given in table [T] and their average speeds are plotted on a 
log-linear scale in figure p] The data from set 1 are given as closed symbols (blue, purple and green) 
in figure [2] the data from set 2 as open symbols (cyan, orange, gold and red). In figure l3b and d, and 
figure [6] the black dots correspond to the mean wavefront speeds of set 1, and the gray ones to the mean 
speeds of set 2. 

In addition to the data shown in figure[2] we also measured the duration of each of the cycles (figurelTk) . 
The numbers we found are consistent with those reported by Foe and Alberts [I] and Parry et al. [20| . 
Averaging over the embryos in each set, we find that the cycle duration can be reasonably approximated 
by a quadratic dependence on the cycle number (figure Wp) . 
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Figure 6. Average speed of each of the two sets of data, on a log-linear plot, fitted by an exponential 
V ^ 2^^("^"°\ e = 0.5 ± 0.05. The black dots correspond to the mean wavefront speeds of set 1, and the 
gray ones to the mean speeds of set 2. 
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Figure 7. Duration of the measured cycles, a) Experimental data. The different symbols and colors 
correspond to the ones in figure pj b) Cycle duration averaged over all experimentally observed embryos 
(black and gray dots for sets 1 and 2 respectively). The cycle durations can be fitted reasonably well by 



a weak exponential i„ = toe 



0.29n 



where to — 33s (set 1) and to = 25s (set 2). 
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D Exponential gradient model 

In this section we investigate the 'exponential gradient timing model' for Drosophila embryo nuclear 
divisions. The idea behind this model is that the nuclei do not signal each other in any way that it is 
time to start mitosis; instead, they sample their local environment for a given protein, such as Bicoid 
(Bed) , and have the length of their cycle depend on the local concentration. Because the concentration 
of many such proteins does indeed exhibit a gradient from one of the two poles, this could explain how 
mitosis starts close to the poles, and then seems to 'travel' along the embryo, where in reality there is no 
traveling front at all. 

Roughly speaking, there are four kinds of patterns expressed by proteins in the early Drosophila 
embryo: an exponential profile along the AP axis starting from one of the poles (of which Bed is an 



example I35|p6]), an exponential profile along the DV axis (such as Dorsal 51 ), a terminal morphogen 
profile which is high at both poles and both low and flat in between (as for the phosphorelation gradient 
of MAPK 1511), and a striped pattern (for e.g. Hunchback, Giant, Paired and Runt [37]). Because 
the observed mitotic waves start at the poles and then spread along the AP axis of the embryo, the 
striped patterns and the DV-axis gradients cannot be the ones causing them. There are often two mitotic 
waves, which start at the two opposite poles, which suggest that the terminal morphogens might be 



good candidates, but they hardly show any profile in the mid-60% of the embryo 51 . The most likely 
candidates are therefore the proteins that have an exponential profile along the AP axis, although in this 
case there must be at least two proteins that can trigger mitosis. This latter observation is a first weak 
point of the model discussed here, but not necessarily cause to rule it out. 

During interphase, the nuclei in the syncytial embryo are surrounded by a nuclear membrane (also 
known as the nuclear envelope) . One of the things the nucleus can do is to concentrate proteins inside 



that membrane. This has been observed directly for Bed by Gregor et al. 36 and confirmed by our own 
observations. Irrespective of whether the proteins are concentrated in the nuclei during interphase, or 
they are distributed throughout the entire cytosol, they always exhibit an exponential decay along the 



anterior-posterior axis. As observed by Gregor et al. 36 , the total amount of Bed steadily increases over 



time, as more of the protein is translated in each cycle. In particular, the amount of protein keeps steady 
pace with the number of nuclei, such that at the start of each cycle, the actual amount of protein in each 



nucleus at a given position in the embryo is always the same 36 . 

Let us denote the position along the long axis by x, and the total length of the axis by L. The local 
concentration at x is then given by c{x) = coe~^'^, where A is the characteristic length scale of the 
exponential profile. As stated above, the experimental results of Gregor et al. tell us that A is the same 



in all cycles, whereas cg goes up 36 . Assuming that all nuclei are equally good at collecting material 
from their environment, the amount of material collected by a single nucleus in a simple one-dimensional 
model of the embryo is given by 

.X+L/2N / T \ 

C{x, N)^ j ^^^ c{y) dy = 2coA sinh (^^^ j e-^/\ (4) 

where N is the number of nuclei. As stated above, the key assumption of the model discussed here is 
that the duration of the cycle of each nucleus depends somehow on the concentration, or rather, on the 
amount of material in the nucleus, so we have 



Aicycic = f{C(x,N)) = f ('2coAsinh f^^ e'-^A 



(5) 



Unfortunately, we do not know what the function / in equation ([5]) is. The only thing we do know 
is that it is monotonously increasing with its argument (the total amount of material in a nucleus) . We 
will therefore explore two explicit possibilities: 

• The simplest possible dependence: / is a linear function. 
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• The dependence that gives the observed behavior of (effective) wavefronts, i.e., that the result- 
ing 'speed' V of mitosis events through the embryo is well-defined and constant throughout, and 

Aicyclc = x/v. 

For the first option, we write f{a) — to — ra, with tg some offset time and r a timescale. We can 
then calculate the speed of an observed mitotic wavefront, as a function of the number of nuclei N, by 
calculating the time difference between two positions x and y in the embryo: 

f AT\ y^'-^ 



Aicycic(y,A^)-Aicycic(a;,iV) 
y-x 



_ (Q) 

2coATsinh(2^) (e-^/^-e-a/A)- ^^ 

Equation mi shows that v depends on the position, so there is no well-defined wavespeed in this model. 
This is of course no big surprise - we just took a random functional dependence for /, so there should 
be no reason to expect it would produce a wavespeed that is position-independent. However, this does 
illustrate the point that a constant wavespeed is something special: we need to specifically choose f such 
that a constant speed comes out. 

Note that it may of course be that there is no constant wavespeed, but that it only appears to be 
constant within our error bars. Although we can not rule this option out, this also does not come 
out naturally. For instance, inserting numbers for Drosophila from Gregor et al. |37| (L = 450/im, 
A — 70/xm — 0.15L, N — 50) we find that for x — Q, the measured speed more than doubles as we take 
the measuring point y across the embryo, which can certainly not be confused with a constant speed. 

To get a constant speed, we need to choose a different function /, specifically a logarithm: /(a) = 
to + T log a. In this case we find: 

In this case, we do indeed find a constant value of v across the embryo. However, we also find that v is 



independent of N. v does depend on the decay length A, but the value of A does not change 36 . The 
wavefront speed predicted by this model is therefore the same for all cycles, in direct contradiction to 
what we observe. 

E Diffusion model 

The process of diffusion is governed by the diffusion equation, here given for a concentration field c(af, t): 

DV'^c{x,t), (8) 



dc{x, t) 2 



dt 

where D is the diffusion constant. Equation ([8| is linear, so we can use the superposition principle: 
the sum of any two solutions is itself a solution. The general solution for a system with no boundaries 
depends only on the initial condition c(x, 0), and is given by 

c(f,t)=y"G(.f,y,i)c(y,0)dy, (9) 

where G{x,y,t) is the Green's function of the diffusion equation, which for a system in n dimensions is 
given by 

1 / |f-y|2^ 



G(x,y,t) = - -— exp -' ^ . (10) 

The Green's function describes the concentration field at x at time t due to a single delta-function 
concentration source at y at time 0. 
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F Mechanical model 

As described in the main text above, we can describe the medium in which the nuclei live as an overdamped 
elastic medium. Motion in this medium can be described by some displacement vector u from a fixed 
reference position. We get force balance by equating the damping forces acting on the nuclei (due to 
friction with the cortical actin layer surrounding the yolk or the outer membrane, and drag due to the 
viscous fluids the nuclei and their surrounding microtubule baskets are immersed in) to the elastic forces 
in the polymer cytoskeleton: 



Tdtu,, = 



E 



2(1 + i^) 



djdji 



E 



2(1-^) 



didjt 



(11) 



As also pointed out above, equation ( |11[ ) is reminiscent of the diffusion equation: a time derivative 
on the left equals second-order space derivatives plus a source term on the right, and it comes as no 
surprise that the solution depends on a quantity with the dimension of a diffusion constant D — A'/F, 



where /i is the material's shear modulus and equals E/(2 + 2v). Moreover, equation (11) also allows 



for a Green's function type solution, but here in the form of a tensor Gijk{x,t), relating an arbitrary 



input Qij5{x)Q{t) to a resulting displacement vector Uk{x,t) 32 . In two dimensions, the input tensor 



Qij has three independent components, and can be decomposed in a (hydrostatic) expansion/contraction, 
a torque and a force dipole. The dividing nuclei in the Drosophila system can be well described by the 
force dipole, with the mitotic spindles pulling the nuclei apart generating the forces, while the volume 
remains fixed and there is no rotation. An example displacement field due to a single force dipole at the 
origin is given in figure |8] 
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Figure 8. Example of a displacement field due to a single force dipole located at the origin and having 
an angle of 7r/6 with respect to the x-axis, obtained by taking the t — >■ cx3 limit of Gijk{x, t)Qij. 
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G Wavefronts 

Now that we know the solutions to the chemical and mechanical diffusion equations due to a single 
source, we can exploit the fact that the chemical and stress diffusion equations (Island 11) are linear to 



compute the behavior of a system with many sources using the superposition principle. For simplicity, we 
pre-arrange the nuclei on a triangular grid, with a little noise in the position of each nucleus to prevent 
artifacts due to a perfect arrangement. This is consistent with our observation that just before the mitosis 
waves the nuclei in an actual embryo have a rather high degree of triangular order, except where there 
are defects due to the fact that a nucleus did not divide in an earlier cycle. Alternatively, we can also 
consider packings with short-range correlations but no long-range order (like the packing of soft repulsive 
spheres) , which gives the same results [32; . 

G.l Wavefronts in the biochemical diffusion model 

We will describe the release by a nucleus of a biochemical with a Dirac delta function source located at 
the position of the nucleus. Because integration is a linear operation, adding two sources, even if they 
divide at different times, is trivial - we simply carry out the integration in equation (|9| for each nucleus 
that has already divided with the time properly offset, then sum over these nuclei. The only problem 
is to determine when each nucleus is supposed to divide. To find out, we perform what is essentially a 
numerical integration over time. We start with a release of material at the origin, which we model by 
having a delta function concentration there at i = 0. By construction, the concentration field is then 
given by G{x, 0, t) as long as no other sources have released their chemicals into the system. We proceed 
in small timesteps At, calculating for each timestep the concentration at the location of each of the 
nuclei that have not yet divided, given the total concentration field generated by the nuclei that have 
divided so far. Suppose there is a total number of N nuclei, M of which have already released their 
chemicals. The zth source is located at x — di, released its chemicals at i = U, and the ii's are ordered. 
For tM < t < tjv/+i we then find by using the superposition principle: 

c(f,t) = y i -expf- '^'^'H : (12) 

where we have taken n = 2. From c{x, t) we can determine when the next source will release its chemicals, 
by solving c{aj,t) = a ioi M + 1 < j < N. We thus check all nuclei that have not yet released anything, 
and determine which one will be the next source by finding the one with the smallest value of t, which 
sets Im+i- If there is a nonzero delay time between activation and release, we simply add it to the found 
value of tM+i- Given the positions of the sources, the system thus has three parameters: the diffusion 
constant D, the release concentration a, and the delay time tdeiay The 'release- wave' is then the time 
at which a given source releases its chemicals versus its distance to the origin (i.e., the first source). An 
example is given in figure |3^, which reveals a clear wavefront with a well-defined wave speed. Figure [9] 
shows the effect of adding delay times: after an initial transit due to finite-size effects near the origin, 
there is again a well-defined wave speed, which slows down with increasing idoiay, and also becomes less 
dependent on the source density as tdoiay goes up. 

G.1.1 Algorithm for finding wavefronts 

In summary, we use the following algorithm to numerically find the wavefronts within the biochemical 
diffusion model (and, with the proper adaptations, for the stress diffusion model as well): 

• Generate a grid of hexagonally arranged nuclei with some small positional noise, centered at the 
origin. 
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• Start with a delta function concentration at the origin at t = 0. 

• Increase time in steps of At. For each timestep, calculate the local concentration at each of the 
sites of the nuclei, due to the nuclei that have released chemical so far. If one of these exceeds the 
critical concentration a, add a delta function concentration peak at this location and time. 

• Stop after either all nuclei have divided or a predefined time has been reached. 



G.1.2 Analysis of the steady-state wavefront 

In the case without delay time, it is fairly easy to determine the speed of the wavefront for the regime 
in which the wavefront is well-established, i.e., when its curvature is small. Suppose the lattice spacing 
is a, and the time it takes to get from one row to the next is t. The speed, in a triangular lattice, is then 
V — a * }j\/3/t, because the spacing between the two layers is a * |v3. To find the time, it turns out to 
be sufficient to only consider the 2 nearest neighbors in the previous row. We choose coordinates such 
that the neighbors are located at (±|,0), and our next nucleus is at (0, |-\/3). Then the time at which 
this next nucleus is activated is given by the solution of 



ATTDt 



exp 



ADt 



(13) 



The results obtained using equation ( 13 ) are almost identical to those obtained by numerical solution of 
the full equations. Corrections can of course be made by including additional sources, but in this case 
that is not necessary. 



To analyze equation ( 13 ) further, we introduce dimensionless variables a = Ana'^a and r = Dt/a? 



The equation then becomes 



a 



-1/4t 



(14) 



which can of course not be solved analytically but is easy to solve numerically. We note that the right 
hand side of (14) has a maximum value of cimax = 8/e at t = 1/4, giving the condition that there can 
only be a wave if a < 2/nea? . If we write the inverse of ( 14) as r = /(a), we can write for the wavefront 
speed 

a 2 /(a) 

so w ex D if both a and a are fixed, but unfortunately the scaling of the speed with a and a is hidden in 
f{a). Based on the numerical determination of f{a) we can capture its features fairly accurately with 
the following function 

4"^ 



/(a) w ha 



Ijr. 



a 



(16) 



where fitting gives hi = 0.93, 62 — 0.14, n — 6.6 and m = 3.6 (see figure 11). Our numerical solutions 



of the full equations show that the data do indeed collapse onto the curve described by equations (15) 



and (16 1 for different diffusion coefficients and nuclear spacings (figure 12). In particular, we find that 
the wavefront speed v always increases as the nuclear spacing a decreases, so v increases with increasing 
cycle number. 



G.2 Wavefronts in the biochemical diffusion model with delay time 



As indicated in Section G.l[ we can include a delay time tdciay into our diffusion model. Now a nucleus 
divides (i.e., releases its chemicals) a time tdeiay afte.r the local concentration first reaches the threshold 
value a. For small delay times, this will simply slow down the wavefront speed a little. However, for 
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even moderately larger delay times, this has severe consequences, because the chemicals released earlier 
will continue to travel outward, pushing the next line closer to the threshold value, possibly even over 
it. We will therefore not be able to suffice with as simple a calculation as before, where the only two 
sources that mattered were the nearest neighbors in the previous row. Instead, we need to include many 
previous rows. Fortunately, the nearest-neighbor approximation still holds: a source is to very good 
approximation only influenced by its nearest neighbors in the previous row, plus their nearest neighbors 
in the row before that, and so on. That means that for each source we have a wedge-shape region of 
older sources that can influence it: 2 in the previous row, 3 in the row before, etc. 

We assume we have reached steady-state, and the wavefront is moving with a constant speed v from 
one horizontal row to the next. The distance between the rows is still ^ayS, and the time it takes the 
front to cover this distance is At. As before, we nondimensionalize, defining a = Ana^a, r = DAt/a^, 
and Td — Dtdc\a.y I o^ ■ The local concentration in dimensionless parameters at the location of the next 
source is then given by 



c(t, Trf) = 


= ^ Cn{T,Td) 




n=no 


Cn{T,Td) = 


- ^ ^ 






exp 



•^n^ -f m^ 



4 



4(nT - Td) 



(17) 
(18) 



In equation (18), we sum over all the relevant sources in row n, which were all activated a nondimen- 



sionalized time nr — Td ago. In equation (17) we sum over the relevant rows: starting at ng, the last row 
that was activated, up to either a large enough value rimax (for true steady-state) or the actual value if 
we want to account for finite-size effects. To determine rip, we note that c„ is not defined for r < r^/n, 
therefore uq is given by the smallest value of tiq such that 

a> lim c(t, Td,no) (19) 

TlTd/no 



is satisfied. Figure 13 shows c{T,Td) for two values of r^, showing how different values of r correspond 
to different values of rig. The function c{T,Td) is a smooth function of r, but for large enough values of 
Td, it is not monotonous, unlike the relation between a and r for the case without delay time (figure [ll| ) . 
The inverse of c will therefore not be continuous for Td larger than some critical value, which also means 
that the resulting wavefront speed is no longer a continuous function of the threshold concentration a 



(figure 14) 



G.2.1 Fitting the data with time delay 

Naturally, the more parameters we have, the easier it is to fit any set of experimental data points. In 
the diffusion model with delay time we have four parameters: the diffusion constant D, the grid size a, 
the threshold value a and the delay time ideiay The grid size (i.e. the spacing between the nuclei) is 
measured independently, leaving us with three parameters which we can vary. Reasonable values for 
the diffusion constant for a small chemical in the early Drosophila embryo are D — 5 ~ 100 /im^/s, as 



measured by Gregor et al. 37 . As indicated in the main text above, we cannot fit the experimental trend 
(a wavefront speed that decreases exponentially with cycle number) with fixed values for idoiay and a for 
any value of D within this range (see figure pp) . The only way we can thus fit the experimental data 
within this model is if either (or both) of a and idoiay change with the cycle number. 

We have systematically investigated a number of options, changing a or tdeiay with cycle number. 
Some results are given in table [2J We did not find any result that fit the data in which the numbers 
change in a well-defined way (e.g. the delay time increasing linearly or exponentially with the cycle 
number). Moreover, in the case of variable delay time, we find that in the first cycle (cycle 10), the 
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interaction is between nearest-neighbors as in the model without time delay, but the interaction range 
goes up every cycle, up to 5 rows apart in cycle 13. In the case of variable concentration threshold, we 
need a change of at least an order of magnitude in each cycle to fit the experimental data. Even in 
the case where we allow both variables to change, we keep finding at least one of these two problems. 
Even though we cannot strictly rule out the diffusion model with variable time delay and concentration 
threshold, these results make it very unlikely that this model is actually correct. 

G.3 Wavefronts in the mechanical model 

The analysis leading to wavefront propagation in the mechanical model is described in Ref. [32]. As in 
the case of diffusion, we must set a threshold to determine when a source (a nucleus) is activated. The 
simplest option is to look at the eigenvalues of the stress tensor: if the largest of those (taking absolute 
values) exceeds a certain threshold a, the nucleus is activated. An activated nucleus adds an additional 
force dipole term to the stress field in the system, which in turn of course affects the displacement field, 



as described by the Green tensor solution of equation (111. Note that we assume linear elasticity, so that 
the strain is linear in the displacement, and the stress is linear in the strain. The superposition principle 
therefore holds not just for the displacements but for the strain and stress fields as well. 

We let each nucleus divide (i.e. generate a force dipole) along a random orientation. The implemen- 
tation of the stress-mediated signaling model follows the same pattern as that of the chemical-diffusion- 
mediated signaling model, with the concentration c replaced by the stress tensor aij . An example imple- 
mentation on a grid of 21 x 21 nuclei is shown in figure [15] The figure shows a clear wavefront which has 
a well-defined speed. 

As in the diffusion model, we analyze our mechanical model in terms of dimensionless parameters. 
There is only one quantity in our model that has the dimensions of a speed, namely i^/aT, which means 
that the resulting wavefront speed has to scale linearly with this factor, as indeed it does. We define a 
dimensionless wavefront speed v = ^— u and a dimensionless stress threshold a = a^a/Q, where Q is the 
strength of the force dipole. We can then write 

«=^g(a,i^) (20) 

We determine the function g{a, v) by numerically solving the model. The best numerical fit to the 
experimental data is shown in figure [2]:. As detailed in 32 , we find that it can be well described by the 



following functional form, derived using a similar argument we used to arrive at equation ( 14 ) for the 

diffusion model: 

4(a-l)log(a) 

g{a, v) = . 21) 

1 — V 



Note that the form in equation (21) differs slightly from that in 32 because of the use of the shear 
modulus fi instead of the Young's modulus E in the definition of v. We use equations (20) and (21) to 
fit the experimental data in figure [2}i. 
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Data set 1 


cycle number 


10 


11 


12 


13 


nuclear spacing (/im) 

wavefront speed (/xm/s) 

cycle duration (s) 

mitosis duration (s) 


23.4 ±0.8 
2.9 ±0.9 
600 ± 60 

237 ±8 


18.2 ±0.6 
2.2 ±0.9 
763 ± 97 
231 ±9 


13.2 ±0.3 
1.5 ±0.4 
922 ± 96 
233 ±12 


9.7 ±0.2 

1.0 ±0.2 

1365 ± 100 

240 ± 12 




Data set 2 






cycle number 


11 


12 


13 






nuclear spacing (^m) 

wavefront speed (/im/s) 

cycle duration (s) 

mitosis duration (s) 


18.0 ±1.2 
4.2 ±0.4 
574 ± 139 
197 ±32 


13.5 ±0.3 
2.9 ±0.3 
757 ±150 
194 ± 26 


10.0 ±0.4 
2.0 ±0.4 

1077 ±160 
194 ± 24 





Table 1. Experimental data averaged over the data sets. Data sets 1 and 2 correspond to two different 
sets of measurements, taken on different days. They correspond to respectively the closed and open 
symbols in figure [21 and figure [Tk. 



cycle 


nuclear spacing 


wavefront speed 


tdclay (s) 


a (10-V/im^) 


(idolay,a) 


number 


(/im) 


(/im/s) 


(a fixed) 


(idoiay fixed) 


(s, 10-4//im2) 


10 


23.4 ± 0.8 


2.9 ± 0.9 


1.5 


0.0044 


(5, 0.254) 


11 


18.2 ± 0.6 


2.2 ± 0.9 


5.7 


0.29 


(10, 0.300) 


12 


13.2 ± 0.3 


1.5 ± 0.4 


16.0 


6.9 


(20, 0.685) 


13 


9.7 ± 0.2 


1.0 ± 0.2 


36.1 


31 


(40, 1.057) 



Table 2. Experimentally determined values of the nuclear spacing and wavefront speed (set 1), and 
numerically determined values of the required delay time tdeiay and threshold concentration a to fit the 
experimental data. In column four, £) = 15 /im^/s and a = 5 • 10~^ jixnT^] in column 5, I? = 10 /im^/s 
and idoiaylO s; in column 6, D = 10 fim^ /s and tdoiay is assumed to double in each cycle. None of the 
columns show a systematic dependency of the parameters on the cycle number, making it impossible to 
assign predictive power to the numbers found, or to find a model to explain the dependencies. Note also 
that for the case of fixed delay time (column 5), we need to assume that the threshold value goes up by 
at least an order of magnitude in each cycle. In the case of variable delay time (columns 4 and 6), the 
interactions become very long-ranged in the later cycles, with nuclei up to 5 rows apart triggering each 
other in cycle 13, even though in cycle 10 the interaction only involves nearest neighbors. 
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Figure 9. Wavefronts in the diffusion model with a time delay, a) Distance to center r vs. activation 
time t shown for a system in which there is a delay time between the local concentration reaching the 
threshold value a and the nucleus releasing new material to the bulk. After a characteristic transition 
regime the system once again exhibits a wave front. Inset: Wavefront visualized in two dimensions using 
a color plot. The wave starts in the center (red dot) with a single Dirac delta peak release. The color 
coding indicates when a nucleus releases its chemical to the bulk, going from red through the different 
hues of the rainbow to violet, b) Effect of the delay time on the wavefront speed for three different 
values of the diffusion constant (blue, purple, green). In a system without any time delay, the dynamics 
is dominated by diffusion, and can be collapsed by normalizing the speed with that of the largest grid 
size (top lines). For high delay times, the dynamics become independent of the diffusion and once again 
collapse (bottom lines). In between, the dynamics depend qualitatively on the exact parameters chosen. 
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Figure 10. Plot comparing the results of the diffusion model to the experimental data (black (set 1) 
and gray (set 2) dots). Purple diamonds: diffusion model with fixed diffusion constant of 100 /im^/s 
and no time delay. Blue diamonds: diffusion model with fixed time delay of 6 seconds (top set) and 10 
seconds (bottom set). Red diamonds: best fits for the diffusion model with the time delay scaling 
exponentially with the cycle number. 
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Figure 11. Plot of /(a), defined in equation (15), determined numerically (blue dots) and fitted with 
equation (16) (red line). 
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Figure 12. Dimensionless wavefront speed v = v/{D/a) as a function of the dimensionless 
threshold a — Aira^a in the biochemical signaling model. The figure shows v for several choices of the 
parameters D and a (different symbols). The data all collapse on a single curve, as described by 
equation (15). The dashed line is the fit to that curve of equation (16). 
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Figure 13. Plot of equation (17 1 for (a) Td = 1-0 and (b) Td = 2.0. The value of no decreases with 



increasing t, plotted values are 1 (blue), 2 (red), 3 (yellow) and 4 (green). Note that c is a continuous 
function of r, but not monotonous for larger values of Td- 
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Figure 14. Time r as a function of threshold value a for several values of the delay time r^: 0.25 
(blue), 0.50 (red), 0.75 (yellow), 1.0 (green), 1.50 (purple). Note that for larger values of r^, r has a 
discontinuity because c(r, r^) is no longer monotonous. For smaller values, r is continuous but does 
have a kink. 
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Figure 15. Wavefront from a simulation with 21 x 21 nuclei. The nuclei are arranged on a hexagonal 
grid with random small offsets. The wave starts at the center point which generates a stress dipole of 
unit strength along the x-axis at i = 0. Whenever the absolute value of the largest eigenvalue of the 
stress tensor at another nucleus exceeds the threshold value a, it also divides, adding a unit stress 
dipole in a random direction to the total stress field, a) Distance of the dividing nuclei to the center vs. 
their time, with a linear fit. b) Graphical representation of the 2 dimensional field, with the dipoles 
indicated in the orientation in which they divide, and color-coded according to the time that they 
divide, on a hue scale (red- yellow-green-blue- violet). 



